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ERROR ANALYSIS IN FOURIER METHODS FOR OPTION 

PRICING 

FABIAn CROCCE, JUHO HAPPOLA, JONAS KIESSLING, AND RAUL TEMPONE 


Abstract. We provide a bound for the error committed when using a Fourier method 
to price European options when the underlying follows an exponential Levy dynamic. 
The price of the option is described by a partial integro-differential equation (PIDE). 
Applying a Fourier transformation to the PIDE yields an ordinary differential equation 
(ODE) that can be solved analytically in terms of the characteristic exponent of the 
Levy process. Then, a numerical inverse Fourier transform allows us to obtain the 
option price. We present a bound for the error and use this bound to set the parameters 
for the numerical method. We analyse the properties of the bound and demonstrate the 
minimisation of the bound to select parameters for a numerical Fourier transformation 
method to solve the option price efficiently. 


1. Introduction 


Levy processes form a rich field within mathematical finance. They allow modelling 
of asset prices with possibly discontinuous dynamics. An early and probably the best 


known model involving a Levy process is the Merton (1976) model, which generalises 


the Black and Scholes (1973) model. More recently, we have seen more complex models 


allowing for more general dynamics of the asset price. Examples of such models include 


to 


the Kou (2002) model (see also Dotsis et al. (2007)), the Normal Inverse Gaussian model 


(Barndorff-Nielsen (1997); Rydberg (1997)), the Variance Gamma model (Madan and 
Seneta (1990); Madan et al. ( 1998| )), and the Carr-Geman-Madan-Yor (GGMY) model 
(Garr et al. (2002 2003)). For a good exposition on jump processes in finance we refer 


Cont and Tankov (2004) (also see Raible (2000) and Eberlein (2001)). 


Prices of European options whose underlying asset is driven by the Levy process 


are solutions to partial integro-differential Equations (PIDEs) ( Nualart et al.| (|2001 ); 
Briani et al. (2004); Almendral and Oosterlee (2005); Kiessling and Tempone (|2011 )) 


that generalise the Black-Scholes equation by incorporating a non-local integral term to 
account for the discontinuities in the asset price. This approach has also been extended 


to cases where the option price features path dependence, for instance in Boyarchenko 


and Levendorski (2002) dHalluin et al. (2004) and Lord et al. (2008) 


The Levy -Khintchine formula provides an explicit representation of the characteristic 


function of a Levy process (cf, Tankov (2004)). As a consequence, one can derive an 
exact expression for the Fourier transform of the solution of the relevant PIDE. Using the 
inverse fast Fourier transform (iFFT) method, one may efficiently compute the option 
price for a range of asset prices simultaneously. Furthermore, in the case of European 
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call options, one may use the duality property presented by Dupire (1997) and iFFT to 
efficiently compute option prices for a wide range of strike prices. 

Despite the popularity of Fourier methods for option pricing, few works can be found 
on the error analysis and related parameter selection for these methods. A bound for 
the error not only provides an interval for the precise value of the option, but also 
suggests a method to select the parameters of the numerical method. An important 
work in this direction is the one by Lee (2004) in which several payoff functions are 
considered for a rather general set of models, whose characteristic function is assumed 
to be known. Feng and Linetsky (2008) presents the framework and theoretical approach 
for the error analysis, and establishes polynomial convergence rates for approximations 
of the option prices. For a more contemporary review on the error committed in various 
FT-related methods we refer the reader to Boyarchenko and Levendorskii (2011), that 
extends the classical flat Fourier methods by deforming the integration countours on the 
complex plane, studying discretely monitored barrier options studied in De Innocentis 


and Levendorskii (2014). 


In this work, we present a methodology for studying and bounding the error commit¬ 
ted when using FT methods to compute option prices. We also provide a systematic way 
of choosing the parameters of the numerical method, in a way that minimises the strict 
error bound, thus guaranteeing adherence to a pre-described error tolerance. We focus 
on exponential Levy processes that may be of either diffusive or pure jump. Our con¬ 
tribution is to derive a strict error bound for a Fourier transform method when pricing 
options under risk-neutral Levy dynamics. We derive a simplified bound that separates 
the contributions of the payoff and of the process in an easily processed and extensi¬ 
ble product form that is independent of the asymptotic behaviour of the option price 
at extreme prices and at strike parameters. We also provide a proof for the existence 
of optimal parameters of the numerical computation that minimise the presented error 
bound. When comparing our work with Lee’s work we find that Lee’s work is more 
general than ours in that he studies a wider range of processes, on the other hand, our 
results apply to a larger class of payoffs. On test examples of practical relevance, we 
also find that the bound presented produces comparable or better results than the ones 
previously presented in the literature, with acceptable computational cost. 

The paper is organised in the following sections: In Section we introduce the FIDE 
setting in the context of risk-neutral asset pricing; we show the Fourier representation of 
the relevant FIDE for asset pricing with Levy processes and use that representation for 
derivative pricing. In Section we derive a representation for the numerical error and 
divide it into quadrature and cutoff contributions. We also describe the methodology for 
choosing numerical parameters to obtain minimal error bounds for the FT method. The 
derivation is supported by numerical examples using relevant test cases with both diffu¬ 
sive and pure-jump Levy processes in Section Numerics are followed by conclusions 
in Section m 


2. Fourier method for option pricing 

Consider an asset whose price at time t is modelled by the stochastic process S = (St) 
defined by St = 5'oe^‘, where X = (Xt) G M is assumed to be a Levy process whose 
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jump measure v satisfies 

( 1 ) 


/r\{0} 


min{y^, l}u (dy) < oo 


Assuming the risk-neutral dynamic for St, the price at time t = T — r of a European 
option with payoff G and maturity time T is given by 

n(r,s) = e-^'"E(G(S’T)|5T-r = s) 

where r is the short rate that we assume to be constant and r: 0 < r < T is the time to 
maturity. Extensions to non-constant deterministic short rates are straightforward. 

The infinitesimal generator of a Levy process X is given by (see Applebaum, 2004) 

Eif{Xt+h)\Xt = x)-f{x) 


f(x) = lim 

h-5-O 


h 


( 2 ) 


1 


= if'ix) + / (/(x + y)- fix) - yl|y|<i/'(x)) i^(dy) 


where ( 7 , cr^, u) is the characteristic triple of the Levy process. The risk-neutral assump¬ 
tion on (St) implies 


( 3 ) 


'|y|>i 


e^p(dy) < 00 


and fixes the drift term (see Kiessling and Tempone (2011)) 7 of the Levy process to 

( 4 ) 


1 


7 = r- a — 

2 


[ (e^ - 1 - yl\y\<i) p(dy) 

Jr\{o} 

Thus, the infinitesimal generator of X may be written under the risk-neutral assumption 
as 

( 5 ) 

C^fix) = (r- fix) + ^f\x) + [ ifix + y)- fix) - (e^ - l)/'(x)) i/(dy) 

V 2 y 2 yR\{o} 

Consider g as the reward function in log prices (ie, defined by gix) = G(5oe®)). Now, 
take / to be defined as 


fiT,x) = Eig iXT)\XT-r = x) 

Then / solves the following PIDE: 

(drfiT,x) =C^fiT,x) 

\/(0,x) =gix), (r,x) e [0,r] X M 

Observe that / and IT are related by 

(6) n(r, Soef = e“'’^/(r, x) 

Consider a damped version of / defined by /^(r, x) = e“"*/(r, x); we see that dr fa = 

^-ax qX j 
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There are different conventions for the Fourier transform. Here we consider the oper¬ 
ator T such that 


( 7 ) 


f e“"/(x)dx 

Jr 


defined for functions / for which the previous integral is convergent. We also use / (a;) 
as a shorthand notation of J^[f] (w). To recover the original function /, we define the 
inverse Fourier transform as 

[/] = f 

Jr 

We have that f (x) = f{x). 

Applying T to fa we get = / (‘^ + fa)- Observe also that the Fourier transform 

applied to f{T,x) gives T {—iio) f (t,uj), where 'h (•) is the characteristic exponent of 
the process X, which satisfies E xhe explicit expression for T (•) is 

(8) T (z) = - y) ^ + y - 1 - (e*' - l)z) p(dy) 

From the previous considerations it can be concluded that 

( 9 ) dr fa = ^ (a - ILV) f {u! - ia) 

Now f {oj — ia) = fa (w) so fa satisfies the following ODE 


( 10 ) 


dr fa (t,U!) / • N 

/ f / =W [a - luj) 

fa {r,OJ) 


^ fa ( 0 , Ul) = Qa (oj) 

Solving the previous ODE explicitly, we obtain 
(11) y(r,cu) = e"'^(“—Va(oj) 

Observe that the first factor in the right-hand side in the above equation is E ^ 

(ie, ipi {—ia — uj)), where ifr (•) denotes the characteristic function of the random variable 


( 12 ) 


ipr (uj) = E (tT (iuj)) 


Now, to obtain the value function we employ the inverse Eourier transformation, to 
obtain 


(13) 
or 

(14) 


1 


fa{T,x)=F~ \fa\ {t,x) = — 6 fa{T, u)duJ 


2 r+oo 

fa{r,x) = - Re 

^ Jo 


27r 

— ZiOX 


fa{T,Uj) 


du 


As it is typically not possible to compute the inverse Eourier transform analytically, we 
approximate it by discretising and truncating the integration domain using trapezoidal 
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quadrature (13). Consider the following approximation: 

Aa; 


(15) 


(16) 


/«,A‘*^,n(r,x) = ^ e (r, ( + | ) Aw 


Acj 


TT 


k=—n 

n—1 

ER' 

k=0 


(k + ^^ Aw 


Bounding and consequently minimising the error in the approximation of /(r, x) by 

/Aa;,n(F)^) = 6 fa,Auj,ni'T~jX') 

is the main focus of this paper and will be addressed in the following section. 


Remark 2.1. Although we are mainly concerned with option pricing when the payoff 
function can be damped in order to guarantee regularity in the sense, we note here 
that our main results are naturally extendable to include the Greeks of the option. 
Indeed, we have by © that 

(17) f(t,x) = ^ [ V, (r, w) dw 

Jr 

so the Delta and Gamma of the option equal 

(18) A {t, x) = = ^ [ {a-iuj) (r, w) dw 

ox 2-71 

(19) r {t, x) = ^ X (r, w) dw 

Because the expressions involve partial derivatives with respect to only x, the results in 
this work are applicable for the computation of A and T through a modification of the 
payoff function: 


(20) ga,A (w) =ga (w) (a - fw) 

(21) ga,r (w) =ga (w) (a - iuf 

When the Fourier space payoff function manifests exponential decay, the introduction of 
a coefficient that is polynomial in w does not change the regularity of ^ in a way that 
would significantly change the following analysis. Last, we note that since we do our 
analysis for PIDEs on a mesh of x’s, one may also compute the option values in one go 
and obtain the Greeks with little additional effort using a finite difference approach for 
the derivatives. 


2.1. Evaluation of the method for multiple values of x simultaneously. The 

East Eourier Transform (EET) algorithm provides an efficient way of computing (15) for 
an equidistantly spaced mesh of values for x simultaneously. Examples of works that 


consider this widely extended tool are Lord et al. (2008); Jackson et ah (2008); Hurd 


and Zhou (2010) and Schmelzle (2010). 


Similarly, one may define the Eourier frequency w as the conjugate variable of some 
external parameter on which the payoff depends. Especially, for the practically relevant 
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case of call options, we can denote the log-strike as k and treat a: as a constant and 
write: 

( 22 ) f 

Jr 

Using this convention, the time dependence is given by 


(23) 


fk,a {r, x) 


g(tuj+a+i)x^^ (cj _ i (a -|- 1)) 

(iuj -|- a) (iuj -|- a -|- 1) 


contrasted with the x-space solution 


(24) 


fk,a {r, x) = 


gpoj a+l)k^^ 


(iuj + a) {iuj -|- a -|- 1) 

We note that for call option payoff to be in we demand that a in (23) is positive. 
Omitting the exponential factors that contain the x and k dependence in (23) and (24) 
respectively, we have that one can arrive from (23) to (24) using the mapping a i—>■ —a—1. 
Thanks to this, much of the analysis regarding the x-space transformation generalises in 
a straightforward manner to the fc-space transform. 


3. Error bound 


The aim of this section is to compute a bound of the error when approximating the 


option price /(r, x) by fa,A‘^,n{'''^x), dehned in (15). Considering 


(25) 


fa,AUr,x) = fr, fk + 

^ fcez 




the total error £ can be split into a sum of two terms: the quadrature and truncation 
errors. The former is the error from the approximation of the integral in (13) by the 


inhnite sum in (25), while the latter is due to the truncation of the infinite sum. Using 


triangle inequality, we have 


(26) £ := |/(r,x) - fAui,n{'^^^)\ 

with 

£q = \U{t,x) - /a,Aa;(r,a:^)| 

£j^ = e“^ |/a,Aa;(T,a:) - /«,Ac,;,n(ua:)| 

Observe that each £, £q and £jr depend on three kinds of parameters: 

• Parameters underlying the model and payoff such as volatility and strike price. 
We call these physical parameters. 

• Parameters relating to the numerical scheme such as a and n. 

• Auxiliary parameters that will be introduced in the process of deriving the error 
bound. These parameters do not enter the computation of the option price, but 
they need to be chosen appropriately to have as tight a bound as possible. 

We start by analysing the quadrature error. 
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3.1. Quadrature error. Denote by Aa, with o > 0, the strip of width 2a around the 
real line: 

Aa ^ {z e C: |Im [z] I < a} 

The following theorem presents conditions under which the quadrature error goes to zero 
at a spectral rate as Auj goes to zero. Later in this section, we discuss simpler conditions 
to verify the hypotheses and analyse in more detail the case when the process X is a 
diffusive process or there are “enough small jumps.” 

Theorem 3.1. Assume that for a > 0." 

HI. the eharacteristie function of the random variable Xi has an analytie extension 
to the set 

Aa — ai = {z ^ C: |Im [z] + a\ < a} 

H2. the Fourier transform of ga{x) is analytic in the strip Aa and 
H3. there exists a continuous function 7 G L^(M) such that fa{T,Lv + i/d) < 7(w) 
for all uj gR. and for all j3 G [—a, o] 

Then the quadrature error is bounded by 

Ma,aiT,x) 


Sq < 


where Ma,a{T,x) is given by 

Ma,a{T,x)\= V [ 

Qr- f -I M 

^£{-a,a} 


(27) 


27r (e2™/A“ - 1 ) 


dw 


Ma,a{x, x) equals the Hardy norm (defined in (28}) of the function ui 1 —)• e a;+ 

i/d), which is finite. 


The proof of Theorem 3.1 is an application of Theorem 3.2.1 in Stenger (|1993|), whose 
relevant parts we include for ease of reading. Using the notation in Stenger (1991^, H\^ 
is the family of functions w that are analytic in Aa and such that 


(28) 

where 


Ircllrri := lim / |tc(z)| d |z| < 00 


Aa{e) = <^zGC: |Re[ 2 ;]|<-, |Im [z]| < a (1 — e) 


Lemma 3.2 (Theorem 3.2.1 in Stenger (1993)). Let w G H\^, then define: 
(29) 


(30) 


I (w) = / w (x) dx 

Jr 

N 

J{w,h) = h ^ w {jh) 

j=-N 

C {w, h) = I {w) — J {w, h) 


(31) 
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then 

(32) 


1C {w,h)\ < 


IP 




2 sinh (^) 


Proof of Theorem \3.1\ First observe that HI and H2 imply that the function w{z) = 

^-i.z+Au/2f^ 

(r, 2 + Alo/2) is analytic in Aa- H3 allows us to use dominated convergence 
theorem to prove that ||ip|| rri is finite and coincides with Mq, ^(t, x). Applying Lemma 

_ Aa ’ 

3.2 the proof is completed. □ 


Regarding the hypotheses of Theorem 3.1, the next propositions provide simpler con¬ 
ditions that imply HI and H2 respectively. 

Proposition 3.3. If a, a and v are such that 

(33) f e*'“^“^^u((iy) < oo and f < oo 

Jy>i Jy<-i 

then HI in Theorem \3 . 1\ is fulfilled. 

Proof. Denoting by (pi (•) the characteristic function of Xi, we want to prove that z i—)• 
ipi{z + ai) is analytic in Aa. Considering that (pi{z + ai) = the only non¬ 

trivial part of the proof is to verify that 


(34) 


p{z,y)i 2 {dy) 


is analytic in Aa, where p: x M —)• C is given by 

p{z, y) = -l-{ey -1) {iz - a) 

To prove this fact, we show that we can apply the main result and the only theorem 
in Mattner (2001), which, given a measure space {Q,A,p) and an open subset G C C, 
ensures the analyticity of f f{-,uj)dp{uj), provided that / : G x D —)• C satisfies: f{z, •) 
is A-measurable for all z £ G', f{-,u}) is holomorfic for all lo G Q] and f |/(•,w)| dp{Lj) 
is locally bounded. In our case we consider the measure space to be M with the Borel 
fj-algebra and the Lebesgue measure, G = Aa and f = p. It is clear that p{x, •) is Borel 
measurable and p{-,y) is holomorphic. It remains to verify that 




\p{z,y)\v{dy) 


is locally bounded. To this end, we assume that Re [z] < b (and, since z £ Aa, Im [z] < a) 
and split the integration domain in |y| > 1 and 0 < |?/| < 1 to prove that both integrals 
are uniformly bounded. 

Regarding the integral in \y\ > 1, we observe that 
(35) \piz, y)\ < + l + + Ifia + a + b) 

for y < —1 we have while for y > 1 we have 

Using the previous bounds and the hypotheses together with 0 and ([^, we obtain the 
needed bound. 

For the integral in 0 < |y| < 1, observe that, denoting f{z,y) = \p{z,y)\, we have, 
f{z, 0) = 0 for every z, dyf{z, 0) = 0 for every z, and \dyyf{z, y)| < c for z G Aa, Re [z] < 
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b, \y\ < 1. From these observations we get that the McLaurin polynomial of degree one 
of y I—>• f{z,y) is null for every z, and we can bound f{z,y) by the remainder term, 
which, in our region of interest, is bounded by obtaining 


(36) 


/0<|2y|<l 


\piz,y)\ z^(dy) < - 


/0<|j/|<l 


y^i/(dy) 


which is finite by hypothesis on z^, which finishes the proof. □ 

Proposition 3.4. If for all b < a, the function x i—>■ e^^^^ga(x) is in L^(M) then H2 in 
Theorem 3.1 is fulfilled. 


Proof. The proof is a direct application of Theorem IX. 13 in Reed and Simon (1975) □ 

We now turn our attention to a more restricted class of Levy processes. Namely, 
processes such that either > 0 or there exists A G (0, 2) such that C (A) defined in 
(37) is strictly positive. For this class of processes, we can state our main result explicitly 
in terms of the characteristic triplet. 

Given A G (0, 2), define C (A) as 


(37) 


C (A) = inf,>i 


y‘^v{dy) 


Jo<\y\<l 

Observe that C (A) > 0 and, by our assumptions on the jump measure p, C (A) is finite. 
Furthermore, if A G (0, 2) is such that 

1 r 


(38) 


lim inf 

eiO 


/ y > 0 

^ J0<\y\<e 


then C (A) > 0. To see this, observe that (38) implies the existence of eo such that 

y‘^v{dy) 



> 0 


0 <|j/|<e 


If eo < 1 observe that 


inf 

eo<e<l 


\ [ y^^^(dy) > > / y'^i^idy) > 0 

^ Jo<\y\<e I Jo<\y\<eo 


'0<\y\<e 

where for the first inequality it was taken into account that ^ > 1 and that the integral 
is increasing with e. Combining the two previous infima and considering |k| = - we get 
that G (A) >0. 

Furthermore, we note that for a Levy model with finite jump intensity, such as the 
Black-Scholes and Merton models that satisfy the first of our assumption, C (A) = 0 for 
all A G (0,2). 


Theorem 3.5. Assume that: a and a are such that (33) holds; ga G and either 
> 0 or C (A) > 0 for some A G (0, 2). Then the quadrature error is bounded by 


£q < 


; Ma,a(r, x) 

27r (eWAo^ - I) 
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where 


(39) M,,,(r,x)= e“"e"'^(“)|fc(ca)| / e 

ce{-i,i} 




duj 


Furthermore, if a'^ > 0 we have 


(40) 




M„,,(r,x)<^ V e“"e"'^(“)|fc(ca)| 
o-Jt , 

ce{-i,i} 


Proof. Considering ha,a{T,x,uj) defined by 


(41) 

we have that 


ha,a{T,X,U!) = ^ 6 (t, W + ZCo) 


ce{-i,i} 


Ma,a{T,x)= / ha,a{T,X,U))duj 
JR 


On the other hand, for j3 G (—a, a): 

(42) |e-'("^+'«"/«(r,a; + i/3)| =e^HA(r,oj + i/3) 


= e^^ 


r'^{c^-\-l3—iu)) 


\ga (U) + if. 


For the factor involving the characteristic exponent we have 

("43^ gT^Ffa+Zl-ia;) _ ^TKe[^{a+l3-iu))\ 

Now, observe that 

( 2 \ 2 
r - yJ + y (^(a + (if - 

(44) + [ (e^^+^^^cos{-yoj)-l-ia + f3){ey-l))u{dy) 

4 r\{o} ^ ^ 

If |a;| < 1 we bound cos{—yuj) by 1, getting 

( 2 \ 2 

r - y J + y ((a + pf - 

+ [ fe{»+h)y _ 1 _ (c, + /3) _ 1)) j,(dy) 

4 ir\{ 0 } ^ ^ 

2 

= 4'(a + /3) - ya ;2 


(45) 
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Assume |a;| > 1. Using that for \x\ < 1 it holds that cos(x) < 1 — x^/4, we can bound 
the first term of the integral in the following manner: 


/ 


g(a+/3)y gQg 


E\{0} 


( 2 /a;)p(dy) < / (l - u{dy) 


+ j e 


< 


(«+/3)?/^(dy) 

[ e^^+l^^yu{dy) - 

Jmxfo} 




'0<|y|<l/|a;| 


y'^i'idy) 


(46) 


< 


4e\{0} 4 


Inserting (46) back into (44) we get 


( 2 \ 2 
r - y j + y (*^“ + 

+ f fe(<-+P)y _i_^a + P){ey- 1 )) u{dy) 
4e\10I ^ ^ 


/E\{0} 

L >| 2 -A 

-Lyc(A) 

2 I |2_ X 

Taking the previous considerations and integrating in M with respect to w, we obtain 


( 1 ^. 

Finally, observing that C (A) > 0 and bounding it by 0, the bound (40) is obtained 
by evaluating the integral. □ 


Remark 3.6. In the case of call options, hypothesis H2 implies a dependence between 
the strip-width parameter a and damping parameter a. We have that the damped payoff 
of the call option is in (M) if and only if a > 1 and hence the appropriate choice of 
strip-width parameter is given by0<a<a — 1. A similar argument holds for the case 
of put options, for which the Fourier-transformed damped payoff is identical to the calls 
with the distinction that a < 0. In such case, we require a < —a. 

The case of binary options whose payoff has finite support {G (x) = '^\x-,x+\ (^)) 
can set any a G M (ie, no damping is needed at all and even if such damping is chosen, 
it has no effect on the appropriate choice of a). 

Remark 3.7. The bound we provide for the quadrature error is naturally positive and 
increasing in Aca. It decays to zero at a spectral rate as Aw decreases to 0. 


3.2. Frequency truncation error. The frequency truncation error is given by 


= 


^Aw 


vr 




k=n 




T, ( A; -b ^ 1 Aw 
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If a function c: (a;o,oo) —)■ (0,oo) satisfies 


(47) 


Re 


+ 0 Au;) <c( 


for every natural number k, then we have that 






< 


TT 


TT 


E 

k=n 

oo 

E' 

k=n 


Re 


-i(fc+|) Atja: j 


k +-]Aw 


1 




k + 


Aw^ 




Furthermore, if c is a non-increasing concave integrable function, we get 

/v/'T’ /*rv^ 

(48) £r<— c(w)dw 

^ J nAu} 


When g'a G and either > 0 or C (A) > 0, then the function c in (47) can be 

chosen as 


(49) 


C(w) = ll^allL- 


[IJQ.OO) 


er'F(a)g- 


E,., 2 ^h!-A 


C(A)1,, 


To prove that this function satisfies (47) we can use the same bound we found in the 


proof of Theorem 3.5, with /? = 0, to obtain 


2-A 


-C^(A)1m>i 


2 I I 

Re [T (a — ioj)] < T (a) — ~ 

from where the result is straightforward. 

3.3. Bound for the full error. In this section we summarize the bounds obtained for 
the error under different assumptions and analyse their central properties. 

In general the bound provided in this paper are of the form 

M 


(50) 




\ g27ra/Aa; _ 


+ 


poo 

/ c (w) dw 

J nAuj 


where M is an upper bound of Ma,a {t,x) defined in (27) and c is non-increasing, inte¬ 
grable and satisfies (47). Both M and c may depend on the parameters of the model 
and the artificial parameters, but they are independent of Aw and n. Typically one 
can remove the dependence of some of the parameters, simplifying the expressions but 
obtaining less tight bounds. 

When analysing the behaviour of the bound one can observe that the term correspon¬ 
dent with the quadrature error decreases to zero spectrally when Aw goes to 0. The 
second term goes to zero if nAw diverges, but we are unable to determine the rate of 
convergence without further assumptions. 

Once an expression for the error bound is obtained, the problem of how to choose the 
parameters of the numerical method to minimise the bound arises, assuming a constraint 
on the computational effort one is willing to use. The computational effort of the numer¬ 
ical method depends only on n. For this reason we aim at finding the parameters that 
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minimise the bound for a fixed n. The following result shows that the bound obtained, 
as a function of Au, has a unique local minimum, which is the global minimum. 


Proposition 3.8. Fix a, a, n, and A and consider the bound £ as a function o/Aw. 
There exists an optimal Aw* G [‘^, 00 ) such that £ is decreasing in (‘^,Aw*) and 
increasing in (Aw*, 00 ); thus, a global minimum of £ is attained at Aw*. 

Furthermore, the optimal Aw is either the only point in which Aw 1 —>■ p (nAw, b) — 
c(nAw), with p defined in (51), changes sign, or Aw = ifp{ujQ,b) — c(wo) > 0. 


Proof. Let us simplify the notation by calling y = nAw, b = 27ran and £ = Tre~°‘^£. We 
want to prove the existence of y* '■ y* > wq such that £{y) is decreasing for ujq < y < y* 
and increasing for y > y*. We have 

M /■“ 

The first term is differentiable with respect to y and goes to 0 if y —?■ O'*'. This allows us 
to express it as an integral of its derivative. We can then express £{y) as 


£{y) = T(wo) + 





ob/uj 


- 1 )' 


w^ 



d{io) 


The first term on the right-hand side of the previous equation is constant. Now we 
move on to proving that the integrand is increasing with y and it is positive if y is large 
enough. Denote by 


(51) 


p{y,b) 


bMe^ly 

[e’^/y - l)^y 2 


Taking into account that c is integrable, we can compute the limit of the integrand in 
00 , obtaining 

M 

lim p{y, b) - c{y) = — > 0 
y^+00 0 

Let us prove that p{y, b) is increasing with y for all b > 0, which renders p{y, b) — c{y) 
also increasing with y. The derivative of p with respect to y is given by 


dyp{y, b) 


bMe^/y {{b/y)e^/y - 2e^/y + b/y + 2) 

y3 [fPly _ 1)3 


in which the denominator and the first factor in the numerator are clearly positive. To 
prove that the remainder factor is also positive, observe that xe^ — 2e^ + x + 2 > 0 if 


3.4. Explicit error bounds. In the particular case when either cr^ > 0 or C (A) > 0 
for some A G (0, 2) we can give an explicit version of (50). Substituting M by M defined 
in Theorem |3.5| and c by the function given in (49) we obtain 


( 52 ) 


£ = £q F £f 
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where 

(53) 

(54) 


^Q= 

ce{-i,i} 


£f = — 

TT 


Cgr4/(ca) \g^{^ca)\ 


TT 


27ra 

eAuj — I 


I 

JR 




POO 

/ ' 

J nA(jj 




-C(A)l|„|>i 


doj 


dw 


This reproduces the essential features of Theorem 6.6 in 


bound (54) can be further improved by substituting 


Feng and Linetsky (2008), the 


L 2 ° by ll^alll,” • 

M [nAcj,oo) 


Remark 3.9. Observe that the bound of both the quadrature and the cutoff error is 
given by a product of one factor that depends exclusively on the payoff and another 
factor that depends on the asset dynamic. This property makes it easy to evaluate the 
bound for a specific option under different dynamics of the asset price. In Subsection 


4.4 we analyse the terms that depend on the payoff function for the particular case of 
call options. 


Remark 3.10. From (53) it is evident that the speed of the exponential convergence of 


the trapezoidal rule for analytic functions is dictated by the width of the strip in which 
the function being transformed is analytic. Thus, in the limit of small error tolerances, 
it is desirable to set a as large as possible to obtain optimal rates. However, non- 
asymptotic error tolerances are often practically relevant and in these cases the tradeoff 
between optimal rates and the constant term becomes non-trivial. As an example, 
for the particular case of the Merton model, we have that any finite value of a will do. 
However, this improvement of the rate of spectral convergence is more than compensated 
for by the divergence in the constant term. 


The integrals in (53) and (54) can, in some cases, be computed analytically, or bounded 


from above by a closed form expression. Consider for instance dissipative models with 
finite jump intensity. These models are characterised by cr^ > 0 and C (A) = 0. Thus 
the integrals can be expressed in terms of the cumulative normal distribution 


(55) 

(56) 


i: 


2 dcu = 


■2 doj = 


2tt 


Ta^ 


2tt 


Ta^ 


1 — ((^V 


Now we consider the case of pure-jump processes (ie, = 0) that satisfy the condition 
C (A) > 0 for some A G (0, 2 ). In this case the integrals are expressible in terms of the 
incomplete gamma function 7 . First, let us define the auxiliary integral: 

I (a,b) = + a “67 f 


for a, 6 > 0. Using this, the integrals become: 


(57) 


. 1 ‘^r 


-C(A)i|„|>i =2(1 + 7 


rC(A) 


2 - A 
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(58) 


r 




^, 2 -A)+l-? ?<1 






2 - A 


? > 1 


An example of a process for which the previous analysis works is the CGMY model 
presented in Carr et al. ( 2002[ 2003), for the regime Y > 0. 

Lastly, when both C (A) and are positive, the integrals in (53) and (54) can be 
bounded by a simpler expression. Consider the two following auxiliary bounds for the 
same integral, in which ? > 1: 


(59) 


i: 


-T[ 


_ 2-A 

doj < 


2 poo I 




. /r? 2 -Ac(A) 


= <^e 2 / 


,2-A 


(60) 


/ 


“ -t( g^t^2+M^^c(A) 


= Yc(A) (i _ $(5^?^)) 


We have that 5(?), defined as the minimum of the right hand sides of the two previous 
equations. 


6W = mi„|5e-YA^/(rk:^,2- 




e 2 <I>(Vto2) - 1 + 26(1) 


is a bound for the integral. Bearing this in mind we have 

(61) [ 

Jr 

and 

^-02) J°° + ^ C(A)1|^I>1 

provided that <^ > 1. 


doj <b (?) 


4. Computation and minimization of the bound 

In this section, we present numerical examples on the bound presented in the previous 
section using practical models known from the literature. We gauge the tightness of the 
bound compared to the true error using both dissipative and pure-jump processes. We 
also demonstrate the feasibility of using the expression of the bound as a tool for choosing 
numerical parameters for the Fourier inversion. 
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4.1. Call option in variance gamma model. The variance gamma model provides 
a test case to evaluate the bound in the pure-jump setting. We note that of the two 
numerical examples presented, it is the less regular one in the sense that cr^ = 0 and 


C (A) = 0 for 0 < A < 2, indicating that Theorem 3.5 in particular is not applicable. 
The Levy measure of the VG model is given by: 


lA/G (dy) = dy ly>o 


Ke-'^+y 


- 1 


Ke^-y 


y 


y<o~ 


y 


and the corresponding characteristic function is given by eq. (7) of Madan et al. (1998): 


(Pr(uj}= ( 1 - + 




K = x 


y- = 


-1 


/6»2y2 




-1 


_ / 6»x 


-1 


By Proposition |3.3| we get that 

(63) a < min {ry_ — a, + a} 

which, combined with the requirement that ga £ (K) (cf. Remark 3.6), implies: 

(64) a < min { 77 + — a,ri- + a, a — 1} 

a < min { 77 + — a,r]- + a, —a} 


for calls and puts, respectively. We note that evaluation of the integral in (13) is possible 


also for a G (0,1) and for a < 0. In fact, there is a correspondence between shifts in the 
integration countour and put-call parity. Integrals with a < 0 give rise to put option 
prices instead of calls. For an extended discussion of this, we refer to Lee (2004) or 


Boyarchenko and Levendorskii (2011), in which conformal deformation of the integration 


contour is exploited in order to achieve improved numerical accuracy. 

In Lee| ( [2004 ) and in our calculations the parameters equal 77 + = 39.7840, rj- = 20.2648 
and K = 5.9311. 

Table presents the specihc parameters and compares the bound for the VG model 
with the results obtained by Lee (2004). Based on the table, we note that for the VG 


model presented in Madan et al. (1998) we can achieve comparable or better error bounds 


when compared to the study by Lee. 


To evaluate the bound, we perform the integration of (27) and (48) by relying on 


the Glenshaw-Curtis quadrature method provided in the SciPy package. To supplement 
Table for a wide range of n, we present the magnitude of the bound compared to the 
true error in Figure 

In Figure we see that the choice of numerical parameters for the Fourier inversion 
has a strong influence on the error of the numerical method. One does not in general 
have access to the true solution. Thus, the parameters need to be optimised with respect 
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n 


Figure 1. The true error and the error bound for evaluating at the 
money options for the VG model test case. 


to the bound. Recall that £ = T(a, Aw, o, n) and £ = £{a, Aw, n) denote the true and 
estimated errors, respectively. Keeping the number of quadrature points n fixed, we 
let (ai,Awi,ai) and (a 2 ) AW 2 ) denote the minimisers of the estimated and true errors, 
respectively 

(65) (ai, Awi, oi) = arg inf £ 

(66) ( 02 , AW 2 ) = arg inf £ 

We further let £\ and £2 denote the true error as a function of the parameters minimising 
the estimated and the true error, respectively 

(67) Ti=T(ai,Awi) 

(68) £2 = £ {o'2, AW 2 ) 

In Figure we see that the true error increases by approximately an order of magnitude 
when optimising to the bound instead of to the true error, translating into a two-fold 
difference in the number of quadrature points needed for a given tolerance. The difference 
between £i and the bound is approximately another order of magnitude and necessitates 
another two-fold number of quadrature points compared to the theoretical minimum. 

In Figurej^ we present the true error^for the Fourier method for the two test cases in 
Table We note that while minimising error bounds will produce sub-optimal results, 
the numerical parameters that minimise the bound are a good approximation of the 
true optimal parameters. This, of course, is a consequence of the error bound having 
qualitatively similar behaviour as the true error, especially as one gets further away from 
the true optimal parameters. 


^The reference value to compute the true error was obtained by the numerical methods with n and 
At*; such that the level of accuracy is of the order 10“^°. 
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Figure 2. The true error £ for the two VG test cases presented in Table 
and the bound-minimising configurations (white circle) {a 2 ,£^ 0 J 2 ) for 
the examples. 


Remark 4.1. In practice, the Hardy norm in coefficient M reduces to evaluating an 
norm along the two boundaries of the strip of width 2o. We find that, for practical pur¬ 
poses, the performance of the Clenshaw-Curtis quadrature of the QUADPACK library 
provided by SciPy library is more than adequate, enabling the evaluation of the bound 
in a fraction of a second. 

For example, the evaluations of the bounds in Tabletake only around 0.3 seconds to 
evaluate on Mid 2014 Macbook Pro equipped with a 2.6 GHz Intel Gore i5 processor, this 
without attempting to optimize or parallelize the implementation and while checking for 
input sanity factors such as the evaluation of the characteristic function in a domain 
that is a subset in the permitted strip. 

We believe that through optimizing routines, skipping sanity checks for inputs and 
using a lower-level computation routines this can be optimized even further, guaranteeing 
a fast performance even when numerous evaluations are needed. 


Remark 4.2. Like many other authors, we note the exceptional guaranteed accuracy of 
the FT-method, with only dozens of quadrature points. This is partially a result of the 
regularity of the European option price. Numerous Fourier-based methods have been 
developed for pricing path-dependent options. One might, for the sake of generality 
of implementation, be tempted to use these methods for European options as well, 
correcting for the lack of early exercise opportunities. This certainly can be done, but 
due to the weakened regularity, the required numbers of quadrature points are easily in 
the thousands, even when no rigorous bound for the error is required. 

We raise one point of comparison, the the European option pricing example in Table 
2 of Jackson et al. (2008), which indicates a number of quadrature points for pricing the 
option in the range of thousands. With the method introduced, to guarantee £ ^ 10“^, 
even with no optimisation, n = 64 turns out to be sufficient. 
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K 

80 

90 

100 

110 

120 

12r = 1 

a 

-16.9 

-13.8 

21.6 

29.10 

36.3 

a 

3.33 

6.45 

18.1 

9.77 

3.52 

II 

CO 

to 

^max 

229 

229 

363 

363 

424 

£ 

3.35 X 10-^ 

0.00334 

0.00562 

3.97 X 10-^ 

7.33 X 10-6 


r 

6 X 10“^ 

0.0032 

0.0058 

6 X 10“^ 

1 X 10-^ 

12r = 4 

a 

-13.8 

-13.8 

22.1 

23.7 

29.10 

a 

6.11 

6.11 

17.9 

15.2 

8.75 

II 

00 

^max 

£ 

62.4 

3.99 X 10"“^ 

42.4 

0.00312 

84.9 

0.00398 

126 

3.57 X 10-^ 

126 

1.33 X 10-6 


r 

1.3 X 10-3 

0.0057 

0.0055 

9 X 10-^ 

1 X 10-“^ 


Table 1. The error bound for European call/put options in the VG 
model for select examples. Refererence result £* from 


Lee 


(2004) 


K 

80 

90 

100 

110 

120 

£ 

2.67 X 10-"^ 

3.49 X 10-4 

4.43 X 10-"^ 

5.52 X 10-4 

6.77 X 10-4 

a 

-1.57 

-1.57 

-1.57 

-1.57 

-1.57 

^max 

22.9 

22.8 

22.6 

22.5 

22.4 


0.34 

0.26 

0.21 

0.17 

0.13 

£^ 

6.87 X lO-'^ 

1.90 X 10-3 

2.82 X 10-3 

2.72 X 10-3 

2.29 X 10-3 


Table 2. Numerical performance of the bound for the Kou model, with 
the test case in Toivanen (2007) (see also d’Halluin et al. (2005)) with the 
number of quadrature points set to n = 32. The point of comparison £* 
refers to the correspondi ng bound c omputed with the method described 


in Chapters 6.1 to 6.4 of Lee (2004) 


In the £\ 


the cutoff error has been 


evaluated using a computationally more intensive Clenshaw-Curtis quad¬ 
rature instead of an asymptotic argument with an exponentially decaying 
upper bound for the option price. 


4.2. Call options under Kou dynamics. For contrast with the pure-jump process 
presented above, we also test the performance of the bound for Kou model and present 
relevant results in This model differs from the first example not only by being dis¬ 
sipative but also in regularity, in the sense that the maximal width of the domain Aa 
is, for the case at hand, considerably narrower. The Levy measure in the Kou model is 
given by 


with p + q 


i^Kou (dy) = A {pe 

1. For the characterisation given in Toivanen] (|2007|) the values are set as 


A = 0.1, r = 0.05 r = 0.25 Sq = 100 
p = 0.3445, 771 = 3.0465 772 = 3.0775 
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from the expression of the characteristic exponent (see Kou and Wang (2004)) 

^2 


a 


'if (z) = z (r - — - XC] + 


z^a^ 


+ A 


PVi Qm 


Vi- z V2 + z 


- 1 


it is straightforward to see 

AaC{z€C:lmz€ (-3.0465, 3.0775)} 


This range is considerably narrower than that considered earlier. In the case of trans¬ 
forming the option prices in strike space, the relevant expressions for option prices as 
well as the error bounds contain a factor exponential in k. The practical implication of 
this is that for deep out of the money calls, it is often benehcial to exploit the put-call- 
parity and to compute deep in the money calls. However, in the case at hand, the strip 
width does not permit such luxury. As a consequence, the parameters that minimize 
the bound are near-identical through a wide range of moneyness, suggesting use of FFT 
algorithm to evaluate the option prices at once for a range of strikes. 


4.3. Binary option in the Merton model. For the particular case of Merton model, 
the Levy measure is given by 


PMerton (dy) — 


\/2 


=exp 


vrcT^ 


(.y - rjf 

2 a| 


and the characteristic exponent correspondingly by 


T(z) 


Merton 


= z \ r — 


a^z 




<72^2 


^ 2^2 

+ A ( -l-z 




4-1 


we may employ a fast semi-closed form evaluation of the relevant integrals instead of re¬ 
sorting to quadrature methods. We choose the Merton model as an example of bounding 
the error of the numerical method for such a model. The parameters are adopted from 
the estimated parameters for S&P 500 Index from Andersen and Andreasen (2000): 

So = 100, A = 0.089, cj = 0.1765, r = 0.05, r,- = -0.8898, aj = 0.4505 


In Figure we present the bound and true error for the Merton model to demonstrate 
the bound on another dissipative model. The option presented is a binary option with 
finite support on [95,105]; no damping was needed or used. We note that like in the 
case of the pure-jump module presented above, our bound reproduces the qualitative 
behaviour of the true error. The configuration resulting from optimising the bound is a 
good approximation of the true error. Such behaviour is consistent through the range 
of n of the most practical relevance. 


4.4. Call options. In Subsection 3.4 explicit expressions to bound S are provided. To 
evaluate these bounds it is necessary to compute ||(7 'q,|| ^oo and IlyaII ^oo • According to 
Remark |3.9[ once we compute these values we could use them for any model, provided 
that they satisfy the conditions considered there. 

The payoff of perhaps the most practical relevance is that of a call option. Consider 
g defined by: 


g{x) = {Soe^-K)+ = So[e^-e^) 
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Figure 3. The true error and the bound E for the dissipative Mer¬ 
ton model, for a range of quadrature points n, along with the bound- 
minimising configurations contrasted to the true error. 


for which a selection of a damping parameter a > 0 is necessary to have the damped 
payoff in (M) and to ensure the existence of a Fourier transformation. In this case we 
have 


(69) 

(70) 

and 

(71) 


9a (w) 


5-0 


exp ((1 — a -|- iu>) x) — exp {k + {ioj — a) x) dx 


S'oexp ((1 — a -|- iu) k) 
{1 + ioj — a) {iuj — a) 


\9a{i^)? 


S2^2(l-a)k 


(a^ -|- cj2) ^(1 — a)^ -|- 


It is easy to see that the previous expression decreases as |w| increase. This yields 

o (l-a)fc 

(72) hah- = \9am = 

and 


(73) — |5'«(?)| 

k.oo) 

The maximisation of I^qI in the strip Aa of the complex plane is more subtle. Denoting 
g'ah, p) = gah + ip), we look for critical points that satisfy dri \gh = 0. This gives 

(74) 4h + h (d/oa -|- 2a^ — 2p — 2a + p'^ + l) =0. 

For p fixed, |(/a| has a vanishing derivative with respect to r/ at a maximum of three 
points. Of the three roots of the derivative, only the one characterised by ry = 0 is a 
local maximum, giving us that for call options 

(75) haWLT = \9a{0,p)\ 

“ pe[-a,a] 
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Now, observe that |ffa(0, p)! is a differentiable real function of p, whose derivative is given 
by the following polynomial of second degree: 

(76) p (p) = (p + a — 2pa — — p^) — 2a — 2p + 1 

We conclude that 

(77) UaWi-p = max{|po(0,p)|} 

where B is the set of no more than four elements consisting of a; —o; and the real roots 
of p that fall in (—a, a). 

Remark 4.3. So far, we have assumed the number of quadrature points n to be constant. 
In real life applications, however, this is often not the case. Typically the user might want 
to choose a minimal n that is sufficient to guarantee error that lies within a pre-defined 
error tolerance. 

In such a case, we propose the following, very simplistic scheme for optimising numer¬ 
ical parameters and choosing the appropriate n to satisfy error smaller than e: 

(1) Select n = no and optimize to find the relevant configuration 

(2) See if £q + < e, if not, increase n by choosing it from a pre-determined 

increasing sequence n = Uj and repeat procedure. 

Especially in using EFT algorithms to evaluate the Fourier transforms, we propose 
Uj = 2^no. We further note that typically the optimal configuration for the optimiz¬ 
ing conhguration for nj+i quadrature points does not differ too dramatically from the 
configuration that optimizes bounds for Uj. 


5. Conclusion 


We have presented a decomposition of the error committed in the numerical evaluation 
of the inverse Fourier transform needed in asset pricing for exponential Levy models into 
truncation and quadrature errors. For a wide class of exponential Levy models, we have 
presented an L°°-bound for the error. 


The error bound differs from the earlier work presented in Lee (2004) in the sense that 


it does not rely on the asymptotic behaviour of the option payoff at extreme strikes or 
option prices, allowing pricing a wide variety of non-standard payoff functions such 
as the ones in Suh and Zapatero (2008). The bound, however, does not take into 


account path-dependent options. We argue that the error for the methods that allow 
evaluating american, bermudan, or knockoff options are considerably more cumbersome 
and produce significantly larger errors so that in implementations where performance 
is important, such as calibration, using American option pricing methods for European 
options is not justified. 

The bound also provides a general framework in which the truncation error is eval¬ 
uated using a quadrature method that remains invariant regardless of the asymptotic 
behaviour of the option price function. The structure of the bound allows for a mod¬ 
ular implementation that decomposes the error components arising from the dynamics 
of the system and the payoff into a product form for a large class of models, including 
all dissipative models. On select examples, we also demonstrate the performance that is 
comparable or superior to the relevant points of comparison. 
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We have focused on the minimization of the bound as a proxy for minimizing numerical 
error. Doing this, one obtains, for a given parametrization of a model, a rigorous 
bound for the error committed in solving the European option price. We have shown 
that the bound reproduces the qualitative behaviour of the actual error. This supports 
the argument for selecting numerical parameters in a way that minimizes the bound, 
giving evidence that this selection will, besides guaranteeing numerical precision, be 
close to the actual minimizing configuration that is not often achievable at an acceptable 
computational cost. 

The bound can be used in the primitive setting of establishing a strict error bound 
for the numerical estimation of option prices for a given set of physical and numerical 
parameters or as a part of a numerical scheme, whereby the end user wishes to estimate 
an option price either on a single point or in a domain up to a predetermined error 
tolerance. 

In the future, the error bounds presented can be used in efforts requiring multiple 
evaluations of Fourier transformations. Examples of such applications include multi¬ 
dimensional Fourier transformations, possibly in sparse tensor grids, as well as time¬ 
stepping algorithms for American and Bermudan options. Such applications are sensitive 
towards the error bound being used, as any numerical scheme will be required to run 
multiple times, either in high dimension or for multiple time steps (or both). 
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Appendix A. Truncation bound in Lee’s scheme 
Eor the strike-space transformed option price we have 

f , N ^ Pi (cu-i(a + l)) 

{ioj + a) {ioj -|- a -|- 1) 

In the particular case of the Kou model, we have 
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We bound the first jump term resulting from upward jumps as 

prii prji {rji — a — 1 + iu) 
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from which it follows 
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similar inequality holds for the term resulting from downward jumps giving us that 

r( + T\( PVi(vi-o,-l) ^ „v2iv2+-‘+i) 

|(^i(a;)|<eV V 2 7 2 V (’'2+“+i)^ J 

The term exp can be bounded from above by an exponential term exp {—'juj): 

/ 2 2 \ 

(78) exp (--— ) < exp (— 70 ; + K) 


—a‘^u? + 27 a; — 2K < 0 

the condition for the quadratic to have no solutions is given by 

7 ^ - 2Ka‘^ < 0, 

setting K = ^, 'Y = a giving us the bound in eq. (78). 

This puts us in place to apply theorem 6.1 of Lee (2004), with 
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